Quantum Kinetic Theory of BEC Lattice Gas:Boltzmann Equations from 2PI-CTP Effective Action 
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We continue our earlier work [Ana Maria Rey, B. L. Hu, Esteban Calzetta, Albert Roura and Charles W. 
Clark, Phys. Rev. A 69, 033610 (2004)] on the nonequilibrium dynamics of a Bose Einstein condensate (BEC) 
selectively loaded into every third site of a one-dimensional optical lattice. From the two-particle irreducible 
(2PI) closed-time-path (CTP) effective action for the Bose- Hubbard Hamiltonian, we show how to obtain the 
Kadanoff-Baym equations of quantum kinetic theory. Using the quasiparticle approximation, we show that 
the local equilibrium solutions of these equations reproduce the second- order corrections to the self-energy 
originally derived by Beliaev. This work paves the way for the use of effective action methods in the derivation 
of quantum kinetic theory of many atom systems. 

PACS numbers: 



I. INTRODUCTION 

In many respects, the dynamics of cold atoms in optical lat- 
tices resemble those of electrons in crystals. Cold-atom sys- 
tems exhibit many favorable attributes for studying quantum 
many -body dynamics, such as the absence of defects in the op- 
tical lattice, and the high degree of experimental control over 
all relevant parameters |7, 'z]. In particular, by varying the 
depth of the optical lattice, the superfluid-insulator phase tran- 
sition can be induced. For weakly-confining optical lattices, 
the system has macroscopic quantum coherence, and interest- 
ing matter wave interference phenomena induced by the peri- 
odicity of the lattice have been demonstrated in experiments 
LIjA' J^^- For tightly-confining lattices, the matter-wave co- 
herence is lost, and the system undergoes a transition to the 
Mott-insulator phase f?]. This regime has become also experi- 
mentally accessible |8, 9, 10|. Outside the weakly interacting 
regime, standard mean field techniques are inapplicable to de- 
scribe the evolution of the system, and alternative methods are 
required. 

Motivated by a recent patterned loading experiment 1 1 ], we 
previously adopted a functional effective action approach ca- 
pable of dealing with non equilibrium situations that require a 
treatment beyond mean field theory ( ll ill , hereafter I). We ap- 
plied the CTP functional formalism [Til] and the two-particle 
irreducible (2PI) effective action |13] to the Bose-Hubbard 
Hamiltonian, and derived equations of motion. This method 
enabled us to go beyond the Hartree-Fock-Bogoliubov (HFB) 
approximation I14l ITsl Il6ll and to incorporate nonlinear and 
non-Markovian aspects of quantum dynamics, which underlie 
dissipation and fluctuation phenomena. 

In its pristine form the 2PI-CTP equations of motion for the 
mean field and the two-point correlation function are compli- 
cated nonlocal nonlinear equations, which defy even numer- 
ical solutions for realistic experimental systems with many 
lattice sites. It is obvious that to get more physical insight 
we need ways to simplify this full theory. In this paper we 
continue this investigation with the goal of showing how to 
formulate a quantum kinetic theory iTtI fisil by way of the 



2PI-CTP formalism. For earlier work addressing this prob- 
lem in quantum field theory, see | l^- iz^ . There exist an ex- 
tensive literature on quantum kinetic theory, many addressing 
BEC dynamics with condensate-noncondensate interactions 
p^ , ^i^-i^ . Those relevant to our present discussions are 

Towards this goal, we ask the question when quantum ki- 
netic theory is a reasonable attainable limit of the more com- 
plete theory based on the 2PI-CTP effective action. Physi- 
cally, a kinetic theory regime exists when the system dynam- 
ics has a clear separation of two time (or length) scales, one 
pertaining to the macroscopic scale describing the kinetic mo- 
tion such as the mean free time and the other to the micro- 
scopic scale such as the duration of collision event. Alterna- 
tively, when perturbations induce disturbances of wavelength 
longer than the thermal wavelengths and frequencies much 
lower than characteristic excitation frequencies, standard ki- 
netic theories may give a reasonable description of the sys- 
tem's dynamics. This is the case for weakly interacting gases 
confined by a slowly-varying external potential. For quantum 
systems, when the quantum features of the many -body system 
act effectively only on the microscopic scale (e.g., when one 
can use a quasiparticle type of approximation), quantum ki- 
netic theory can provide an adequate description. It fails when 
such a two-time separation does not exist, such as in strongly 
correlated systems or systems with macroscopic quantum co- 
herence ' . 

The organization of this paper with a brief of our findings 
is as follows. In SeclUwe summarize our prior results for the 



' We have in mind systems whose quantum coherence or correlation or en- 
tanglement extends to macroscopic dimensions. Examples are coherence 
tunneling phenomena 1 4C], quantum properties of microelectro-mechanical 
systems |41L |43| and of course, BEC, which certainly has macroscopic 
quantum coherence. The impossibility of a two-time separation refers only 
to the condensate state alone. The interaction between the condensate and 
the non-condensate atoms can under general conditions allow a two-time 
separation and a kinetic theory description, as is the topic of our present 
discussion and much prior work 
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HFB and second order equations of motion fll'l and express 
them with hghtened notation in a more compact form. In Sec. 
mil we discuss how a quantum kinetic theory can be derived 
from a quantum theory of interacting particles. We first dis- 
cuss this issue under more general conditions, where a two- 
time separation may not exist. A kinetic theory is obtained 
from the full hierarchy of correlation functions by truncation 
of higher order correlations and the imposition of causal fac- 
torizable conditions. We use the «PI-effective action to il- 
lustrate this conceptual framework. In Sec. II VI we focus on 
situations where there is a two-time separation in the system 
dynamics. We delineate the physical conditions and show the 
procedures in deriving quantum kinetic equations from the 
2PI-CTP equations of motion. Then we introduce further sim- 
plifications and discuss how to derive the familiar Boltzmann 
equations. In Sec|V]we study how these kinetic equations ad- 
mit, as a special yet important case, equilibrium solutions. We 
show that under the Popov approximation the second order 
2PI equations yield to the same second-order damping rates 
originally obtained in Beliaev's pioneered work |43] but with 
a modified effective mass due to the presence of the lattice. In 
Sec lVII we conclude with a few general remarks. 



II. 2PI-CTP EFFECTIVE ACTION FOR THE 
BOSE-HUBBARD HAMILTONIAN 

Here we summarize the structure of this method and collect 
the useful equations obtained from our earlier investigation 
lHHl . We will refer to the numbering of equations therein with 
a prefix I. 



A. The Bose-Hubbard Hamiltonian 

The one dimensional Bose-Hubbard Hamiltonian is given 

by 

H ^ -J ^(aj fli+i + a,|+i fti ) + ^ V^a\ai 

i i 

+ i^U^a\a\ata, , (1) 

i 

where hi and a\ (called $i , $| in Paper I) are the bosonic op- 
erators that annihilate and create an atom on the site i. Here, 
the parameter U denotes the strength of the on-site repulsion 
of two atoms on the site i; the parameter Vi (called in Pa- 
per I) denotes the energy offset of each lattice site due to an 
additional external potential that might be present (such as 
a magnetic trap), and J/fi denotes the hopping rate between 
adjacent sites. Next-to-nearest neighbor tunneling matrix el- 
ements are typically two orders of magnitude smaller than 
the nearest-neighbor ones, and to a good approximation they 
can be neglected. The Bose-Hubbard Hamiltonian should 
be an appropriate model|44] when the loading process pro- 
duces atoms in the lowest vibrational state of each well, with 



a chemical potential smaller than the energy separation to the 
first vibrationally excited state. 

The classical action associated with the Bose-Hubbard 
Hamiltonian Q, is given in terms of the complex fields 
and a* by 



S[a*,ai\ — j dt'^^iha*dtai + J dt J (a* 

i i 



U 



To compactify our notation we introduce aj(6 = 1,2) defined 
by ai = aj, a* — af. In contrast to I, where we set Vi (t) — 
0, here we allow the presence of an external potential Vi in 
S[a*, ai]. In the derivation of Boltzmann equations we will 
assume that Vi is a slowly-varying function in position and 
time and treat it as an external perturbation. 

In terms of these fields the classical action takes the form 



S[a] ^ j dtY,\haba1{t)hdta\{t) 

i 

dtY, (jcTaba1+^{t)a\{t) - \v{t)aa,bat{t)a\{t) 

i ^ 

/ dtY,(^-^iaabant)a1it)) 



(2) 



(3) 



where Af is the number of fields, which is two in this case, and 
summation over repeated field indices a,b — (1, 2) is implied. 
hab and aab are matrices defined as 



hab = i 



-1 

1 



O'ab = 



1 

1 



(4) 



In terms of the familiar Pauli matrices, aab = cTx and hab = 
-ay. 

We define the following index lowering convention 

Xa — CFabX^. (5) 

After second quantization the fields a° are promoted to op- 
erators. We denote the mean field or the expectation value of 
the field operator by zf{t) (called 4>f{t) in Paper I) and the 
expectation value of the fluctuation field (fi, ipf = a" — zf by 
Gfj{t,t'). Physically, |zf (i)P is the condensate population 
and the two point functions Gij{t,t') determines the quan- 
tum fluctuations around the mean field: 



G^^(i,t') 



(6) 



The brackets denote taking the expectation value with re- 
spect to the density matrix and Tq denotes time ordering along 
a contour C in the complex plane. 

In order to describe the non-equilibrium dynamics we spec- 
ify the contour of integration to be the Schwinger-Keldysh 
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contour fl^] along the real-time axis or closed time path 
(CTP) contour. Using the CTP contour, the two-point func- 
tions are decomposed as 

G^^lt, t') = e,tp{t, t')Gi';> {t, t') + 9,tp{t', t)Gp {t, t'), 

(7) 

(8) 
(9) 



where 



with Octp {t — t') being the CTP complex contour ordered theta 
function defined in Eq. (1.30). 

All correlation functions of the quantum theory can be ob- 
tained from the two particle irreducible (2PI) effective action 
r[z, G]. In Ref. Ill] we showed T[z, G] is given by: 



T[z,G] 



-TrlnG-' + -TrD-\z)G 
2 2 ^ ^ 



S[z] 

+r2[z,G]+ const, (10) 
where iD~'^{z) is the classical inverse propagator given by 



SS[z] 



5zf{t)5z%t') 



(11) 



= {Sijhhabdf + JiSi+ij + Si^ij)aab) S{t - t') 
{2z,a{t)z,b{t) + aabzt{t)z,c{t)) 5^j5{t ~ t'), 

and r2[z,G] consists of all two-particle irreducible vacuum 
graphs in the theory (the diagrams that do not become discon- 
nected by cutting two propagator lines) with propagators set 
equal to G and vertices determined by the interaction terms 
in S[z + (fi] . 

The dynamical equations of motion for the mean field zf (t) 
and the propagators Gf'{t, t') are found by solving the equa- 



tions = and 5g-(m') 

24) and (1.26) respectively. 



sr[z,G] 



0. They were given in (I. 



The action F including the full diagrammatic series for r2 
gives the full dynamics. It is of course not feasible to ob- 
tain an exact expression for r2 in a closed form. Various ap- 
proximations for the full 2PI effective action can be obtained 
by truncating the diagrammatic expansion forr2. The ones 
relevant for this paper are the HFB approximation and the 
full second order approximation. The HFB approximation 
corresponds to a truncation of r2 retaining only the first or- 
der diagram in U which is z independent, i.e. keeping only 
the double- bubble diagram (Fig. 1 in paper I). The full sec- 
ond order approximation corresponds to a truncation retaining 
also diagrams of second order in [/ ( the basket-ball and the 
setting-sun). 

Hereafter, to lighten the notation, we introduce a more com- 
pact set of symbols for the physical quantities than was used 
in Paper I, which contains more details: 



{am \ ( ziU) 

G<\{tX) = G>;{t\t) 



(12) 

(13) 
(14) 

(15) 

(16) 



The notation ti means that the function must be evaluated at 
the time t and at the lattice site i. 



B. The HFB and full second order equations of motion 

The equations of motion derived from the 2PI-CTP ef- 
fective action( (1-24) and (1-26)) have terms which can be 
grouped as the single particle, the HFB and the second order 
contributions, as follows: 



^ / dt"{D-\u,tl)-S"^^{U,tl))H(i!it'^) dt"S{U,tl)H{tlt'^), 
k •' k '' 

^ j dt"{D-\u,tl)~^"^^{U,tl))g[tlt'^) =Y,j dt"nU,tl)g(tlt'^)-5,,5c{t-t'), 



(17) 
(18) 



r 



where ^ {ti , ) is the inverse free particle propagator given az is the Pauli matrix: 
by: 

D-\t,,t'j) EE {i5,j<j,dt + J{5^+ij+5,^ij))5{t-t') 



-S,,V,S{t-t'), 



(19) 



1 

-1 



(20) 
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gHFB ^jjj YjHFB jjjg HFB self-energies of H and g re- second order corrections. 

spectively, and 5 and S are the remaining parts of the self- Using Eq. (1-43) it can be shown that S^^^ and Y.^^^ are 

energies of H and g, which we will assume are given by the given by: 



^HFB 



{ti, t: 



U 



ij^ {Tr {H{U,t'^) + g{U,t'^)) I + 2 {H{U,t'^) + g{U,t'^))) 5{t - 
. U 



S^'^^iU, t;.) = tj^ {Tr {H{U,t'^) + g[U,t'^)) I + 2g{U, t'^)) S{t ~ t')5,, 

I 



(21) 
(22) 



where / is the identity matrix. 

In paper I, we used the CTP contour of integration (which 
is also usually called "in-in" contour) to evaluate the second 
order contribution. Use of the CTP formalism was impor- 
tant there, because it provided the technical means to formu- 
late our initial value problem in a completely causal manner, 
removing the Feynman boundary conditions on the Green's 
function used in the conventional "in-out" formalism 1 17]. In 
this work we are more interested in deriving kinetic equations 



J 



which are especially devised to study relaxation of systems 
close to equilibrium. With this purpose in mind, as done by 
Kadanoff and Baym iflTiflsil . it is better to set the initial con- 
ditions in the far past. We follow them hereafter and use the 
CTP contour, but instead of setting the initial time to zero, as 
was done in paper I, we choose it to be — cxo. The equations 
of motion we obtain in this way agree with the equations of 
motion and are given by: 



dt" {D-\u,4) - s"^^{u, ti)) H{ti t'^) - / dt"7(t„ ti)H{ti <;.) 

^ -oo J — OO 

/OO pi' 
„ -OO <J —oo 



= 0, 



(23) 
(24) 



/OO 
dt" {D-\t,,tl) - g'^^Kt'lA) - 

Y^f dt"m,t';)g(^\tlt'^)-Y, f dt"E(^\t,,t'DAitlt'^] 

1^ J -oo J-oo 

/OO 
dt"g^^\t,,tl) {D-\tlt'^) - = 

Y^f dt"A{t,,tl)i:(^\t'it'^)-Yf dt"g^^\t,,tl)T{tlt',, 

1, -J — OO 1^ J — oo 



(25) 



(26) 



In the above equations, Eq. \26\ is the hermitian conjugate of per I we denoted these by a subscript ^''^ . Here, for ease of 
Eq. M5\ . Eq. MA\ is the hermitian conjugate of Eq. J23> and comparisons with the literatur^we have changed to the nota- 
we have introduced the spectral functions 



7(i„t;.) = {s>{t,,tr)-s<{t,,t'^)), 
r(t„t;') ^ (s]>(t,,t;')-s<(i„i;')), 

A{t,,t';) EE ig>{t,,t'^)~g<it,,t'l)). 



tion of Kadanoff and Baym llTll . 7, F and A. We will show 
later that 7 and F contain information about the condensate 
(27) and noncondensate particle decay rates respectively. 

(28) 

(29) 



Notice that A{ti, t"), r{ti, t"), j{ti, t") are just the spectral If we use the full second order expansion, (1.60) and (1.61), 
functions defined in (I. 36) multiplied by a minus sign. In Pa- S^^-* and E^^^ are given by 
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+2g(^HU,t'^) (^H{t'^,U)g'^^HU,t'j) + g^^Ht'j,U)H{U,tr) + g^^\t'^,U)g^^^ (31) 
grr [H{t'^,U)g(^\u,g + g^$\t'^,U)H{U,g + g^$\t'^,U)g^^\u, } . 

I 



It is convenient to decompose the above equations in their 
matrix components. To do that we introduce the definitions 



a>(t t') - -if .32^ _.orj2 



-iSU 



At equal times, the quantities and are related by the 5< (^^^ ^ mj-i(m,j-m*, + 2pjip,j), (37) 

bosonic commutation relations. Using Eqs. ( I32t and ( I33t into -'^ 
the self-energy equations we get 



yHFB(. ,f^ _ 2U f 2\z.,\^ + p.u + pa zf + mu \ , 



V 2\Z^\' + p^^ + p, 

2U_ / |zip + pu + pu mu 
My m*^ l^iP + Pu + Pii 



t^) = ^2 (Pj^Pfj + 2mijPijmj* + 2p.ijmj*ZiZj + pfjZjZi* + 2pjiPijZiZj* + (40) 

'^TTlijTTlji Zj Zi ~\- '2inij p-ij Z{ Zj 

^i2{ti:'^j) = {2pjrmijpij + 2pjipijZiZj + uiij'^mji* + 2mijZiZjmji* + 2mijPijZjZi'' + (41) 

2pjiTn^j Z^Zj -\- TJlij Zi Zj ), 

= ^2 (P^j'^Pji + 2pij7njimij* + 2pijZiZjmij* + pj^ZjZ.* + 2pijPjiZiZj* + (42) 
2TTijiZ{inij Zj H~ 2pijTnj{Z{ Zj ), 

^f2(^«'^j) = {2p^jmjiPji + 2pijPj^ZiZj 4- mji^m^j^ 4- 2mjiZ^Zjmij^ + 2p.,jmjiZjZ^'' + (43) 

2nzjiPjiZiZj -\- TTiji Zi Zj ), 



^2^2 (t, , t;. ) = (i, , t;. ) {p,, ^ p,, } , (44) 

5^?2 ! *j ) = (^i ; ) {^i ^ , Pii ^ Pii } , (46) 

S|l(^i7ij ) = i;^2(^i>^i) {^i ^ ^i*>"^j* ^ -(47) 



The above expressions for the self-energy, which contain two- 
particle irreducible diagrams up to second order in the inter- 
action strength, agree exactly with those used in Refs. i26l 
In Ref. |30|] the authors used these equa- 
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tions as the starting point of a quantum kinetic theory be- 
fore applying the Markovian approximation. It is important 
to mention that in contrast to other self-energy approxima- 
tions that may lead to equations of motion that do not sat- 
isfy conservation laws, the 2PI effective action formalism is 
a "^-derivable" i45l l46ll approximation and therefore all the 
equations of motion derived from it are guaranteed to be con- 
serving. Moreover, as we showed in paper I, a truncation up 
to second order in the interaction strength is not appropriate 
to describe far-from-equilibrium dynamics outside the weak 
coupling regime. Away from the weak coupling regime, the 
2PI effective action can be a powerful tool. For example a 
1 /TV expansion of the 2PI effective action has been shown to 
provide a practicable controlled nonperturbative description 
of far-from equilibrium dynamics without the small coupling 
restriction L25J . 



m. FROM QUANTUM THEORY OF INTERACTING 
PARTICLES TO QUANTUM KINETIC THEORY 

From previous sections it can be observed that the equations 
of motion obtained from the 2PI effective action are quite in- 
volved: nonlinear and nonlocal integro-differential equations, 
not readily solvable in closed form. To progress further we 
need to introduce approximations based on physical consider- 
ations. It is easier to proceed if one can observe and justify 
a separation of time scales in the relevant physical processes 
in question, i.e., one related to quantum processes which are 
usually microscopic in scale (note quantum entanglement and 
correlation of the system can extend to much greater scales, 
meso or even macro) and one related to the kinetic or trans- 
port properties, which is usually macroscopic in scale. How- 
ever, this assumption of a scale separation, may not be valid 
in mesoscopic processes (as in strongly correlated systems) or 
macroscopic quantum coherence effects (see footnote 1). For 
those situations where a separation of macroscopic and mi- 
croscopic time scales which would permit an effective kinetic 
theory description does not exist, one can adopt the effectively 
open system framework quantified by the nPI-CTP effective 
action and the hierarchy of equations it generates. We begin 
with a discussion of the latter situation which is more demand- 
ing and general. We describe the conceptual pathway for the 
construction of quantum kinetic theory from the nPI effective 
action. Though somewhat theoretical and formally oriented, 
it may be of some use, as this is the first point of contact with 
quantum kinetic theory from the effective action approach, in 
the atomic and molecular physics (AMO) context(For more 
details on this subject see yj|20|,|21|,|22|], where our discus- 
sions in the following section are based on). 



A. Quantum kinetic theory from («PI) effective action 



(relativistic) Wigner functionfTT*] and its higher-correlation 
analogs, which are obtained by a Fourier transform of the 
relative coordinates in the Schwinger-Dyson equations ll48ll 
for the correlation functions, or alternatively, in the mas- 
ter effective action (defined as the nPI effective action when 
n oo, we are dealing with n = 2 here) whose varia- 
tion yields the Schwinger-Dyson equations. This is a quan- 
tum analogue of the BBGKY hierarchy |49], expressed in a 
representation convenient for distinguishing between micro- 
scopic (quantum field-theoretic) and macroscopic (transport 
and relaxation) phenomena. As such, it does not require near- 
equilibrium conditions, and in fact, is applicable for a rather 
general moment expansion of the initial density matrix 

To understand how quantum kinetic theory is derived from 
an nPI effective action and how it relates to the familiar Boltz- 
mann's theory, it is perhaps helpful to examine the relation 
between this theory in its full generality and an effective 
Boltzmann description of relaxation phenomena for the one- 
particle distribution function of quasiparticles. In nonequilib- 
rium statistical mechanics, as is well known lEollsill . the act 
of truncating the classical BBGKY hierarchy does not in itself 
lead to irreversibility and an i7-theorem. One must further 
perform a type of coarse graining of the truncated, coupled 
equations for n-particle distribution functions. For example, 
if one truncates the hierarchy to include only the one-particle 
and two-particle distribution functions, it is the subsequent as- 
sumption that the two-particle distribution function at some 
initial time factorizes in terms of a product of single-particle 
distribution functions (which is at the heart of the molecular 
chaos hypothesis where the colliding particles are initially in- 
dependent, but correlated after a collision ) what leads to the 
(irreversible) Boltzmann equation. The assumption that the 
two-particle distribution function factorizes is an example of 
a type of coarse graining called slaving of the two-particle 
distribution function to the single-particle distribution func- 
tion, in the language of i20ll . The situation in quantum ki- 
netic field theory is completely analogous. One may choose 
to work with a truncation of the hierarchy of the Wigner func- 
tion and its higher correlation analogs, or one may instead 
perform a slaving of, for example, the Wigner-transformed 
four-point function, which leads (within the context of per- 
turbation theory) directly to the (relativistic) Boltzmann equa- 
tion 1 19] and the usual i7-theorem |22]. Typically this slav- 
ing of the higher correlation function(s) involves imposing 
causal boundary conditions to obtain a particular solution for 
the higher correlation function(s) in terms of the lower order 
correlation functions fl9ll20ll . The truncation and subsequent 
slaving of the hierarchy within quantum kinetic field theory 
can be carried out at any desired order, as dictated by the initial 
conditions and relevant interactions. As with any coarse grain- 
ing procedure, in implementing the slaving of a higher cor- 
relation/distribution function to lower correlation/distribution 



It may be useful to begin by defining what we mean by a 
quantum kinetic theory. It contains, but supersedes, the quan- 
tum version of Boltzmann's theory. Formally it refers to the 
theory based on the hierarchy of coupled equations for the 



^ It should be pointed out that in order to identify the Wigner function with 
a distribution function for quasiparticles, one must show that the density 
matrix has decohered, and this is neither guaranteed nor required by the 
existence of a separation of macroscopic and microscopic time scales l50ll . 
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functions, one is going over from a closed system to an ef- 
fectively open system, the hallmarks of which are the emer- 
gence of dissipation fioll and noise/fluctuations |20]. This 
fact has led some to search for stochastic generalizations of 
the Boltzmann equation |52], motivated by the fact that sys- 
tems in thermal equilibrium always manifest fluctuations, as 
embodied in the fluctuation-dissipation relation. (A deriva- 
tion of the stochastic Boltzmann equation from quantum field 
theory can be found in iSoll .') 

The essential point about the process of slaving of higher 
correlation (or distribution) functions is that it is a step which 
is independent of the assumption of macroscopic and micro- 
scopic time scales. In fact, a completely analogous proce- 
dure exists at the level of the Schwinger-Dyson equations 
(i.e., without Wigner transformation) for correlation functions 
in an interacting quantum field theory |20]. Recall that the 
Schwinger-Dyson equations are, in the context of nonequilib- 
rium field theory in the Schwinger-Keldysh or closed-time- 
path (CTP) formulation, an infinite chain of coupled dynam- 
ical equations for all order correlation functions of the quan- 
tum field. The importance of the closed-time-path formalism 
in nonequilibrium situations is that it ensures that the equa- 
tions are causal and that the correlation functions are "in-in" 
expectation values in the appropriate initial quantum state or 
density matrix. As with the BBGKY hierarchy in nonequilib- 
rium statistical mechanics, the common strategy is to truncate 
the hierarchy of correlation functions at some finite order A 
general procedure has been presented for obtaining coupled 
equations for correlation functions at any order I in the cor- 
relation hierarchy, which involves a truncation of the master 
effective action at a finite order in the loop expansion i20ll . 
By working with an / loop-order truncation of the master ef- 
fective action, one obtains a closed, time-reversal invariant set 
of coupled equations for the first / + 1 correlation functions, 
z ~ Ci, g — C2, C3, . . . , C;+i. In general, the equation of 
motion for the highest order correlation function will be linear, 
and thus can be formally solved using Green's function meth- 
ods. The existence of a unique solution depends on supplying 
causal boundary conditions. When the resulting solution for 
the highest correlation function is then back-substituted into 
the evolution equations for the other lower-order correlation 
functions, the resulting dynamics becomes non- time-reversal 
invariant, and generically dissipative. As with the slaving of 
the higher-order Wigner-transformed correlation functions in 
quantum kinetic field theory, we have then gone over from a 
closed system (the truncated equations for correlation func- 
tions) to an effectively open system. In addition to dissipa- 
tion, one expects that an effectively open system will manifest 
noise/fluctuations (an example of slaving the four-point func- 
tion to the two-point function in the symmetry-unbroken \^'^ 
field theory is given in |20]). Thus a framework exists for ex- 
ploring irreversibility and fluctuations within the context of a 
unitarily evolving quantum field theory, using the truncation 
and slaving of the correlation hierarchy. 

While it is certainly not the only coarse-graining scheme 
which could be applied to an interacting quantum field, the 
slaving of higher correlation functions to lower-order correla- 
tion functions within a particular truncation of the correlation 



hierarchy, as a particular coarse graining method, has several 
important benefits. First, it can be implemented in a truly non- 
perturbative fashion, where the variance of the mean field can 
be on the order of the "classical" mass (defined as the sec- 
ond order derivative of the effective potential in the equation 
of motion for the mean field, which provides the natural time 
scale of the system dynamics). This necessitates a nonpertur- 
bative resummation of daisy graphs (the leading contributions 
in a large N expansion) 1 13], which can be incorporated in 
the truncation/slaving of the correlation hierarchy in a natu- 
ral way. ^ Second, the truncation of the correlation hierarchy 
accords with our intuition that the degrees of freedom readily 
accessible to physical measurements are often limited to the 
mean field and two-point function. 



IV. SYSTEMS WHOSE DYNAMICS ADMIT TWO-TIME 
SEPARATION 

An alternative (actually more common and easier) route to 
reach a kinetic theory description from n-body quantum dy- 
namics becomes available when there is a clear separation of 
two time scales in the system dynamics. This is the usual text 
book treatment of kinetic theory we are familiar with. The two 
different scales in the system are the time (or length) scale sep- 
aration between the duration of a collision event (or scattering 
length) and the inverse collision rate (or the mean free path). 
For quantum processes, in the weakly interacting regime, we 
expect there is also a separation between the kinetic scale of n 
particles (expressed in the center of mass coordinate) and the 
quantum scale (expressed in the momentum corresponding to 
the Fourier transform of the relative coordinates between two 
particles), which describes how quantum processes (such as 
radiative corrections) change the particles' mass-energy and 
momenta. Using these approximations it is possible to recast 
the full quantum dynamics into the simpler forms of two cou- 
pled equations which constitute quantum kinetic theory, the 
Boltzmann equation governing the distribution functions and 
what is known as the gap equation for the modified dispersion 
relation. 

For a three dimensional uniform Bose gas the duration of 
a collision event tq is given by the time that a particle with 
average velocity v spends in the interaction region measured 
by the range of the two-particle interaction potential. This 
range for a repulsive potential is typically given by the s-wave 
scattering length and thus tq ~ a-s/v. On the other hand, 
the inverse collision rate Tc or time between successive col- 
lisions is approximately given by Tc ~ (na^v)'^, where n 
is the particle density. The required separation of time scales, 
Tc 3> To implies the inequality na^ ^ 1 or in other words, the 



' At late times in the thermalization stage, when the quantum field is near 
equilibrium, an effective kinetic description may be justified, but will likely 
require resummation of hard thermal loops (see, e.g., 15311 ^ . Under such 
circumstances, even the evaluation of transport coefficients is nontrivial for 
high temperatures I23tl54ll . 
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necessary condition required for the validity of a scale separa- 
tion is that the system must be in the dilute weakly interacting 
regime. For atoms in optical lattices the dilute weakly inter- 
acting conditions required for the scale separation is fulfilled 
if the average repulsive interaction energy U n, where n is the 
mean number of particles per lattice site, is much smaller than 
J, the quantum kinetic energy needed to correlate two atoms 
at adjacent lattice sites, or C/n/ J <C 1. 

Perhaps an intuitive way to understand the scale separa- 
tion is the following. At equilibrium the correlation func- 
tions describing a homogeneous system are translationaly in- 
variant and stationary. If the system is disturbed from equi- 
librium, collisions among particles would break both invari- 
ances. However, as long as the the interaction energy per 
particle is smaller than the typical kinetic energy per parti- 
cle, inter-particle collisions are few and far between. In this 
case the quantum-mechanical entanglement between collision 
partners decays faster than the time required for the next colli- 
sion to take place, particles can be considered as free between 



where we have introduced the definition A6'(tj+i/2) = 
d{ti^i) — 9{ti). As shown in Ref. y^, in the context of 
the discrete Bose-Hubbard model, it is convenient to map the 
unitary gauge transformation to the so called phase-twist of 
the Hamiltonian. The twisted Hamiltonian exhibits additional 
phase factors, e^'^^ in the hopping term, which are known as 
the Peierls phase factors I55ll56ll . 

The scale separation is performed by introducing the vari- 
ables: 

i?=(z + i)/2, T=(t + t')/2, (51) 



r=ii^j), r={t-t'), (52) 

For a translationally invariant system at equilibrium, the con- 
densate density no{ti) is position and time independent and 
the propagators g'-^'> {ti,tj) only depend on the relative coor- 
dinates variables r and t and are highly peaked about their 
zeros. If the system is disturbed by small perturbations, such 
as an external potential V{ti) which varies slowly in space 
and time, we expect for systems with scale separation, that the 
gauge-transformed propagators, g'-^^ti, t'^), acquire a slowly 
varying dependence on the center of mass coordinates R and 
T but still to be peaked around the zeros of r and t. We em- 
phasis that the gauge transformed, not the original variables, 
are the ones that are expected to be slowly varying. The rea- 
son is that even if the perturbation is slowly varying, the phase 
d{ti ) can be a rapidly varying function and it can induce strong 
variations in the condensate amplitude and in the propagators. 



collisions and approximate time and space translational invari- 
ance holds. 



A. Coarse-graining procedure 

To make the scale separation, for BEC systems at hand, it 
is best to perform first a gauge transformation which makes it 
easier to identify (and coarse-grain away) the fast variations 
induced by the rapid change of the condensate phase. Follow- 
ing Ref. 11811 we introduce the gauge transformation 

= e^'^'^^^UoiU), (48) 



where ynj^ and 6{ti) are real. The equations of motion 
are invariant under the phase transformation if we replace 



^e{u+i/2)g^^^^ + e-*"^'^''(*'-i/^)(5,;_i,)) 5{t - t'), (50) 
I 

Before going further, it is important to discuss the issue that 
by defining the spatial center of mass coordinates at points 
that strictly speaking are not lattice sites points we might be 
introducing un-physical degrees of freedom. We stress though 
that this is not the case for system with scale separation. Un- 
der the slowly varying approximation the un-physical degrees 
of freedom are excluded, since the functions evaluated at the 
R points may be thought of as the average over neighboring 
physical lattice sites. 

We proceed now to describe the coarse-graining procedure 
that uses the slowly varying property of the propagators in the 
center of mass variables to simplify the equations of motion. 

If the phase twist applied to the system is small ^ vr, 
the Peierls phase factors can be written as, e'"^^ = 1 — i/S.9 — 
\^0^. In this case, the phase factors can be physically con- 
nected to the imposition of an acceleration on the lattice and 
the energy change resulting from the phase twist can be at- 
tributed to the kinetic energy of the superflow generated by 
the acceleration. Under this picture in the context of the Bose- 
Hubbard model the quantity can be also connected, as is 
the gradient of the phase in non lattice systems, to the super- 
fluid velocity: 

fiVs{U+r/2) = 2JM{t,+^,2)ai. (53) 

with ai the lattice spacing. 

If the disturbances introduced by the perturbation are small, 
the superfluid velocity is expected to be a slowly varying func- 
tion in space and time and to a good approximation its second 
order variations can be ignored, i.e. Aus(t) = '2,[vs{ti+i/2) — 
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Vs{ti)],y i- Again, the quantity Vs{ti)) may be thought of pendence of the superfluid velocity, the propagator D^^ can 
as the average over neighboring lattice sites: [^5 (^^-1-1/2) + thus be written in terms of the superfluid velocity as: 
t's(tj-i/2)]/2. Using the small angle and slowly varying de- 

I 

D-\U, t'j) « (^S,,[iha,dt - hdt6{U) - V{U) - JvliU)] + J(l + ^a,A«,(t))[(5,+i, + 5{t - t') 

+ {J\a,Vs{Um+ij~5,-i,))5{t-t'). (54) 

I 



where we have introduced the dimensionless superfluid veloc- 
ity ^;,(tO = ^1^. 

At equilibrium, the time derivative of the phase is related 
to the chemical potential. Extending this identification to the 
nonequihbrium system we define the chemical potential as 

ix{U) = -hdtOiU) - Jvl{U) -V{U), (55) 



no{ti) 
V{U) 

Similar approximations can be made on the two point 
functions by introducing a change of variables {ti ,tj) 
(r, r; R, T). The slowly varying dependence in R and T al- 
low us to treat g^^\r, r; R, T) as a continuous functions in 
the center of mass coordinates and neglect second order varia- 
tions in them. On the other hand, it is important to include the 
discrete dependence on the r — i — j variables, inherent to the 
tight binding Hamiltonian, in order to retain all the quantum 
effects introduced by the lattice which are crucial to a proper 
description of the system. 

We now introduce a Fourier transform with respect to the 
relative coordinate variables. Since hereafter we use the 
gauge-transformed functions exclusively, the primes will be 
dropped to simplify the notation: 

g<-^\U,t'^)^g<^^\rT;RT) (60) 

q 

H{U,t))^H{rT-RT) (61) 
= -^^lY. J da;e(^^"'™)i?(i?,(7;T,c.), 



If we make a change of variables (ti) {R + (r/2),T + 
(t/2)) in the one-point functions no{ti), ^{ti), Vs{ti) and 
V{ti) and use the the slowly varying dependence of the func- 
tions on the center of mass coordinates, to a good approxima- 
tion the functions can be treated as continuous functions and 
second order variations in R and T can be neglected. Thus, 
they can be written as: 

(56) 
(57) 
(58) 
(59) 



using Eq. ( 15 6> in Eq. (I61> we get: 



H{R, q; T, Lu) = 27rM (/ + <j,) n,{R, T)S{u;)Sgo. (62) 



In Ea.( l62t . the quantity no{R,T) is just related to the con- 
densate density of atoms at the space time point {Rai,T). In 
Ea. (l60> . the upper diagonal component of the two-point func- 
tion gf^{R,q]T,Lo) corresponds to the well known Wigner 
distribution function |47]. It can be interpreted as the den- 
sity of noncondensed particles with quasimomentum q and 
energy hui at the position Rai and time T. On the other hand, 
g{'i{qR;T,u!) is essentially the density of states available to 
a particle that is added to the system at {Rai , T) with quasi- 
momentum q and energy huj. As opposed to a normal system, 
the presence of the condensate gives nonzero values to the 
off-diagonal terms of the functions gj^^ (i?, T, u>). We refer 
to them as the anomalous contributions to the respective two 
point functions. 



= no{R + (r/2), T + (t/2)) = ^^(i?, T) + ^9fln„(i?, T) + ^-drnoiR. T), 
= ii{R + (r/2), T + it/2)) = ii{R, T) + ^-dRti{R, T) + ^-drt^iR, T), 
= v,{R+ (r/2), T + it/2)) - v,{R, T) + "^dnv^iR, T) + ^drv^iR, T), 
= ViR + (r/2), T + it/2)) = y(i?, T) + ^duViR, T) + ^-drViR, T). 
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B. Generalized Boltzmann equations 

The generalized Boltzmann equations can be obtained as 
the Fourier transform of the equations of motion for the case 
in which the variations in R and T are very small: in par- 



ticular when the inverse propagator ^ and the self energies approximated by; 



vary very little as Rai is changed by a characteristic excitation 
wavelength or T is changed by an inverse excitation energy. 

If we neglect the second order variation in T and R, as ex- 
plained above, the equations of motion Eq.(l23> to ( I26> can be 



-^S+'-^)h = [D-\H] + [^S, H] + \ [7, H] 



H D 



-'-[H,D-'] + '-[H,^S]~\[H,j]. 



(63) 
(64) 



D-\g^^-> 



J 4 L . 



g(^\^j: 



(65) 



(66) 



with 



D-\qR;T,uj) = (a, {hu -v,{R,T)2J ^Yx^^qai)) + (2J cos{qa{) ^ ^liR^T)) I) . (67) 

I 



In Eqs. ( I63l66> all the quantities depend on (q_R; TuS). 
In the equations we have also introduced the following func- 
tions: 



dependence of the quantities on R and T . The brackets in 
Eqs. ( I63I66> denote the generalized Poisson brackets defined 
as: 



^S^{R,q;T,Lj), 



^S{R,q;T,u;) 

sR5](i?,g;T,w) = Y."'' {R, q;T,uj) 
^Y.^{R,q;T,uj), 



^S^{R,q;T,uj) 
^Y.^{R,q;T,uj) 
^giR,q;T,Lu)^P J 



duj' "f{R,q;T,uj) 

27r LU — LU' 

duo'T{qR-T,uo') 
2tt lo — lo' 
dw' A{R,q]T,uj') 
27r LO — lo' 



(68) 
(69) 

(70) 
(71) 
(72) 



with P denoting the Cauchy principal value and 7(i?, q; T, lo), 
T{R,q;T,Lo), S"^ {R,q:T,uj), T."^ {R,q:T,uj) and 
A{R,q;T,uj) understood as Fourier transforms of the 
functions -f{U,t'^), T{U,t'^), s"'' {U.t'j), j:"^{U,t'^) and 
A{ti,t'j) respectively. 

To approximate the discretized equations by the continuous 
differential equations we have also used the slowly varying 



dAdB_ _ dAdB_ 



dnAd.B-d.AdRB. (73) 



Notice that even though the continuous limit has been taken 
at the kinetic scale, the discreteness introduced by the lattice, 
crucial for a correct description of the physics, is taken into 
account at the quantum scale, as can be seen in Eq. ( I67l i 
where the free propagator has a trigonometric dependence on 
the quasimomentum q, characteristic of lattice-type systems. 
If the disturbances in the system are small enough that only 
long wavelength modes are excited, qa ^ 1, the excitations 
only see the lower quarter of the band. In this case the free 
propagator reduces to 



D-\{q « l/ai)R; T, io) ^a,ihio~ Vs{R, T)p) + 
2J ~ ^ + H{R,T)] I . (74) 



2m* 



which is like the free propagator for a non-lattice system and 
the role of the lattice is just to introduce an effective mass m* . 
Here p — hq and m* — Y? j {2afj). 
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If we define the statistical functions(which carry a super- 
script in Paper I, Eq. (1.35)) as: 



F{R,q;T,uj) 
U{R,q;T,co) 



g>{R,q;T,u;) + g<{R,q;T,Lo) 
2 

^>iR,q;T,Lo) + i:<iR,q;T,u) 



(75) 



(76) 



Eq. ( I65> and Eq. ( I66t . can be rewritten in terms of statistical 
and spectral functions as: 



F D-^ - SRS - -r 



^9 - 2^ ) n 



D-'-m + -r,F 



^g-'-A,Il 



(77) 
(78) 



- KE) A - mg = -- { - A] - [r, Kg]} , 



^-{[A,D-'-^^~[^g,T]}. 



r 



(79) 
(80) 



Eqs. ( I63t . (I64> and ( I77> - ( I80> are our passage to the Boltz- 
mann equations. They describe the state of the gas at a given 
time. Different from the HFB equations they include colli- 
sional integrals for binary interactions. 



effects due to interactions, i.e. how interaction effects change 
the energy momentum dispersion relations from that of free 
particles to a more complicated spectrum. Because these two 
effects are physically distinct, we can treat the left and the 
right hand sides in a different way. 



C. Ordinary Boltzmann equations 

To progress further we can introduce more simplifications 
based on physical considerations. The ordinary Boltzmann 
equation emerges from the approximation in which the self 
energies that appear on the left side of Eas.(l63>. i64\ and (I77> - 
( I80> are handled differently from those which appear on the 
right. These two appearance of the self-energy play a different 
physical role in the description of the dynamics IITII . The self 
energies on the right hand side describe the dynamical effects 
of collisions, i.e., how the collisions transfer particles from 
one energy-momenta configuration to another On the other 
hand, the self energies on the left describe the quantum kinetic 



In the derivation of the ordinary Boltzmann equations, one 
completely neglects all the kinetic effects in the second or- 
der self energies (the dependence on T and R in the second 
order self-energy terms on the right hand side) and retain dy- 
namical effects (T and R dependence on the left hand side). 
In this way, we get the familiar Boltzmann equations which 
describe the particles as free particles in between collisions 
with a modified energy-momentum dispersion relation. It is 
a reasonable assumption in dilute weakly interacting gases in 
which the duration of a collision is very short compared to the 
essentially interaction-free dynamics between isolated colli- 
sions. Neglecting kinetic effects in the second order self ener- 
gies, Eas.( l63> . (I64> and (I77> - ( I80> can be approximated to 



D-i _ + l^j H = [D~^ S"^^,H] , (81) 

H ^D-' _ - ^7) = [H, - S"^^] , (82) 

D-^ - + Ir) F - n {^g + \a^^-\ [^o ' - ^"^"^ .F] , (83) 

-m:)A-mg = -'- [D-^ - a] , (84) 

F [d-^ - - Ir) - {^^g - 1a) n = -1 [F, - E^^^] , (85) 

A ~ - sRgr = - ^ [A, D'^ - E^^^] . (86) 
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If we take the trace of the sum and the difference of each 
one of the above equations with its hermitian conjugate, they 
can be simplified to : 



Moreover, if we define the operator IrM = M12 



M. 



21 



and apply it again to the sum and the difference of each one 
of the equations Ml\ to ( I80> with its transpose we also get: 



Tr{{D-^ - F - nSR.g} 
Tr {{D-^ - ^J:) A - mg} 
Tr[D-^ ^S"^^,H] 



0, 
0, 
0, 



(87) 
(88) 
(89) 



-Tri-fH), (90) 
--Tr{TF -UA)X9l) 
0. (92) 



Re {Ir {{D-^ ~ ^S) H}) 
Re (> {{D-^ - 3fiS) F - n5R.g}) 
Re {lr{{D-^ - KE) A - mg}) 

Im {Ir {{D-^ - nS) H}) 
Im {lr{{D-^ - 3fiS) F - li^g}) 

im {^^r{{D-^ -m:)A- mg]) 



hm (> [D-^ - S^^^, F] + >(rF - UA)) , 
\lm{-\r [D-'-i:"^^,A]), 
-\Re{-\r [D-'-S"^^,H]+^{^H)), 
-^Re {Ir [D-^ - S^^^, F] + lr{TF - UA)) 
-\Re{lr [D-'-J:"^^,A]), 



(93) 
(94) 
(95) 
(96) 
(97) 
(98) 



r 



with Re and Im denoting the real and imaginary parts. To 
close the set of equations, we need an equation of motion for 
the superfluid velocity which can be found from the defini- 
tions Ea. (l55> and Ea.( l53> to be: 



Eas. ( 1871981 together with Eq.(l99> form a closed set of 
equations that describe the state of the gas at a given time. 
Equations ( I87I89I I and ( I93I95> are usually called the gap 
equations. They describe the quantum properties of a gas 
which is evolving according to Boltmaznn-type equations (|90j- 
ED and ^6l%\ . Under the derived formalism Eqs.(|H7)-(|55) 
form a coupled set of equations which replace the original 
dynamics. The equations have to be solved self consistently 
for any analysis. 



EQUILIBRIUM PROPERTIES FOR A HOMOGENEOUS 
SYSTEM 



There are two situations in which we expect an equilibrium 
solution to come from the Boltzmann equations. Firstly when 
the system has never been disturbed and it remains in its equi- 
librium state. Secondly when the system has had sufficient 



time to relax after an applied perturbation. In this section we 
will show how the second order nonequilibrium Boltzmann 
Equations lead, in these special cases to the linear equilib- 
rium solutions obtained from the HFB approximation fl6ll up- 
graded with second order corrections in U . 

At equilibrium, in the absence of any external potential, the 
functions and H are completely independent of R and T. 
In this case the generalized Poisson-bracket terms are zero and 
Eqs (|89l|92j and (|95l|92} imply that: 



A{D-^ ~ KE) - (Kg)r 



0. 



(100) 



Because 5R5(g, uo) is determined by A{q^ lu) as indicated in Eq. 
M2\ . Eq. ( I100> is satisfied when A{q, lo) is given by 



iA{q,uj) 



- + -r 



D-^ - - -r 

2 



(101) 



and the function 'Stg{q, uj) given by 
duj' A{q,uj') 



UJ — UJ' 

■ 1 

- KE + 



D-^ - - -r 



(102) 
-1 ~1 
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From Eqs.j90ll, ( I96> . ( I91> and (|^} we also get, at equilibrium, 
the conditions 



TF 



7 
HA 



(103) 
(104) 



Eqs. ( I103> and (I104> are just the mathematical statement of 
detailed balance. They represent the physical condition that at 
equilibrium the net rate of change of the density of particles 
with momentum q and energy uj is zero. Since it is always 
possible to write iflTII 



F{q,uj) - [n,{Lu) + ^]Aiq,u), (105) 



Eq.( ll04> can only be satisfied if 



n(g,L^) 



rigiuj) 



r{q,to) 



(106) 



is satisfied. Detailed study of the structure of the self-energy 
indicates that nq{uj) is related to the Bose-Einstein thermal 
distribution, = ^i3^_^ with (3 interpreted as the local 

inverse temperature in energy units fT7l,fl8'|. In refs. i'st'i the 
authors prove that the only translational invariant solution is 
the thermal. 

Since H contains delta functions in momentum and energy at 
equilibrium, we get from Eq. (I87> : 

= -2J + KS'ii(0,0) + 5R5i2(0,0). (107) 



A. Quasiparticle formalism 

In the noninteracting case the diagonal terms of A{q, cu) are 
just delta functions with peaks at values of hu that match the 
possible energy difference which results from adding a single 
particle with quasimomentum q to the system. In the many 
body system the energy spectrum is sufficiently complex so 
that the diagonal elements of A{q, lu) are not delta functions 
but instead continuous functions of uj. However, there are al- 
ways sharp peaks in A. These sharp peaks represent the co- 
herent and long lived excitations which behave like weakly 
interacting particles. These excitations are called quasiparti- 
cles. From Eq. jl01> it is possible to see that the quasiparticle 
decay rate is determined by F. The quasiparticle approxima- 
tion is obtained by considering F very small for small values 
of oj. This assumption implies that D^^ = D^^ — — |F 
is essentially real with only an infinitesimal imaginary part. 
The zeros of about which A is very sharply peaked are 
identified with the quasiparticle energies fkUq. 

Using the assumption of a very small F, and the identity 

lim =P — ^ -inSiuj-uj'), (108) 

e^O U! — UJ + le U! — Uj' 

it is possible to write the matrix components of D^^ as: 



D- 

huj 



fl 

-1 



(109) 



Cqq{q,Uj) Mq-q{q,Uj) 

M*{-q,-Uj) £.* J-q,-uj] 



with 



Cqq^uj) — — 2 Jcosqa; — /i 



^HFB 

^11 



duj' Tii{q,Lu') 



2tt uj — uj' 



Mq-q{uj) = Sf/^(g,C^) + 



duj' Ti2{q,uj') 
27r UJ — uj' + ie 



(110) 



(111) 



The quasiparticle amplitudes Uq and Vq are the solutions to the 
eigenvalue problem 



^qqiq,^q) Mq-q{q,UJq) 
M*q^q{-q,-UJq) C*qq{-q, -Uq) 



Uq 
-Vq 

(112) 
P = 1. 



and satisfy the normalization condition \uq\'^ - 
In the absence of vortices it is always possible to find an 
ensemble in which the amplitudes {uq,Vq) are purely real 
and Uq = U-q, Vq — V-q. In terms of the quasiparticle 
amplitudes, the matrix elements of the spectral function A, 
Eg. ( lion , are given by: 



Aiiiq,uj) 



-2Im 



UJ — UJq + iO+ 



UJ — UJq + iQ 



Ai2{q,uj) = 21m. 



2^[ul5{uj-ujq)-vl5{uj + ujq)\, (113) 

v„u„ 



'■q <Jq 



uqu,q 



UJ — UJa 



^0^ 



uj„ + iO 



A22{q,Uj) 

A2i{q,i^) 



— 2nUqVq[d{u; — UJq) — S{UJ + UJq)] , (114) 

-Aiii-q,-u;), (115) 
Al2{q,u). (116) 



Finally, using the definitions of F and A, we can express the 
matrix components: Pq{uj), Pq{uj) and mq{uj) defined as the 
Fourier transform of p,y, /Jy and m,y respectively (see Eqs. 
\2>2\ and \2>2>V ) in terms of quasiparticle amplitudes: 



Pqiuj) = 2TT[u^^nq{uj)5{uj - UJq) 

Wq{l + nq{u))5{uj + Uq)] , 
Pq{uj) = 2Tl[u^q{l nq){uj)5{uj - UJq) 

mq{uj) 



+Vqnq{uj)S{uj + UJq)] , 
2TTUqVq [nq{u})6{uj + U}q) 
-{1 + nq{uj))S{uj + UJq)] 



(117) 

(118) 
(119) 



B. HFB approximation 

Under the HFB approximation the matrix SftE and SftS" are 
just given by Y.^^^ and S^^^ . In terms of the quasiparticle 
amplitudes and setting Af = 2, they can be written as: 
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gHFB 



U 
U 



2 (rio + n) no + ni 
Uo + fh 2 (rio + n) 

rio + 2n m 
m Tio + 2n 



(120) 
(121) 



with 



m 



(2^5(^5) + 1). 



(122) 
(123) 



In the HFB approximation. Eg. (II 12> and Eq. J107> then yield: 



— 2 Jcos(qa;) — /i + 2U{no + n) 
U {uo + fh) 



U {uo + ni) 
—2Jcos{qai) — ^ + 2U{no + h) 



hujq { ^ 



(124) 



-2J+[/no + 2C/n + C/m. 



(125) 



As a final step, to fix the total density to n, the constraint 

n = Uo + h, (126) 

has to be satisfied. 

For a given density and temperature Eqs. J124t - J126t form 
a closed set of equations. At zero temperature, they reduce to 
the HFB equations derived in 1 16| using the quadratic approx- 
imation. 

The Hugenholtz-Pines theorem state s is^ that a homoge- 
neous system at equilibrium has to fulfill 



£,,(0,0)-X,_,(0,0) = 



(127) 



The above equation implies that the energy spectrum of a Rose 
gas is gapless,i.e. there is an excitation with an energy that 
tends to zero in the limit of zero momentum. Mathematically 
the theorem implies that the two-point propagator g{q, ui) has 
a pole at q = u! = 0. Physically it reflects the fact that small 
rotations of the phase of the condensate wave function cost 
little energy (Goltstone mode of the broken symmetry). The 
Hugenholtz-Pines theorem is a consequence of the invariance 
of the mean field and the two point propagators under a phase 
transformation. 



The HFB approximation violates the Hugenholtz-Pines the- 
orem: 



£,,(0, 0) - Mq-q{0, 0) = -2Um ^ 



(128) 



One way to solve the gap problem is to set the anomalous 
term fh to zero in HFB equations. This procedure is known as 
HFB-Popov approximation. The HFB-Popov equations were 
first introduced by Popov |59], and at equilibrium they are 
consider a better approximation than the HFB equations be- 
cause they yield a gapless spectrum. Nevertheless the HFB- 
Popov equations are not conserving and therefore they are not 
appropriate to describe dynamical evolution. 



C. Second-order and Beliaev approximations 



When second order terms are taken into account the matri- 
ces Cqq and A4q-q become energy dependent. For simplicity 
we restrict the calculations to the zero temperature case when 
Uq = 0. In terms of the quasiparticle amplitudes the contribu- 
tions to the self-energy at second order are given by 



15 



hM 



2[/2 ^ 



k,p 

Cqq{q,Uj) = 

hM 



= Utlo + Um + 

/ 2AkBq-k + 2CkAp + 2CkBq-k + 5CkCq-k 2BkAq-k + 2CkAp + 2CkBq-k + SCkCq-k 

\ U! -UJk~ i^q-k + ie UJ + UJk+ LUq-k " K 

2AkBpCq—k—p ~^ ^k^p^q—k—p '^BkApCq—k—p ^" ^k^p^q—k- 



(129) 



■E 
E 



UJ - UJk - UJp - l^q-k-p + «e + LJfe + CJp + iOq-k-p " 

-2Jcosqai - fi + 2Uno + 2Un+ (130) 
A A-fc + 2Agg-fc + 4Cfc^,_fc + 2CkCq-k BkBq-k + 2(gfc A-fc) + 4Cfcg,-fc + 2CfcC,_fc ' 



H = -2J + 



AkApBq—k—p ^ 2^/1/^CpCg— /e— p 



UJ + ujk+ ujq^k — ie 

BkBpAq—k—p ^" 2BkCpCq — k—p 



C/?io + 2Un + Um - — — V , 

k,p 



le J \ + cj/j + cjp + Uq^k-p — «e 
/ 2^fcSpCfc+p + 2BkApCk+p + 2CkCpCk+p 



2lP 



E 

k,p 



UJk +UJp+ UJk+p 
2AkCpCk+p ^" AkApBk^p ~\' 2BkCpCk-^p BkBpAk-\-p 



UJk+ UJp + LOk+p 



where the quantities A, B and C are defined as 



Ak = Ufe, Bk =vl, Ck = -UkVk 



(131) 



The inclusion of second order terms modifies the structure 
of the HFB equations. The matrix that we need to diago- 
nalize to find the quasiparticle energies depends now on the 
quasiparticle mode in consideration. This means that a sepa- 
rate nonlinear problem must be solved for every quasiparticle 
state, whereas the solution of the HFB equations yields the 
whole quasiparticle spectrum. The matrix which is to be di- 
agonalized also becomes intrinsically nonlocal and to solve 
for a quasiparticle state with quasimomentum q we have to 
sum over alldifferent quasimomenta. Finally, the diagonal el- 
ements are no longer equal as was always the case in all the 



quadratic approximations. 

If we omit the second order terms containing no condensate 
amplitudes, the equations that we get are the tight-binding ver- 
sion of the ones originally derived by Beliaev |43]: 



with 



and 



Mq-q{q,Uj) = Uno + XAMq-qiq,Uj) 



^qq{q,^) = 



A* = 



eM + 
-2 J 



Urio 
-Uur, 



\ACqq{q,Uj) 

f AA/i. 



e[q] = AJs\T?{qai/2) 
A^i = 2Un + Um. 



(132) 
(133) 
(134) 



(135) 
(136) 



2U^ 



>E 



AMq^q{q,uj) = Um + —no ^ , 

fe ^ 

2BkAq-k + 2CkAp + 2CkBq-k + iCkCq-k 



( 2AkBq^k + 2CkAp + 2CkBq-k + '^CkC, 



q—k 



UJ - UJk - ujq^k + ie 



ACqq{q,uj) = -Um 



uj + ujk+ ujq-k - «e 

2U^no f AkAq-k + 2AkBq-k + ^CkAq-k + 2CkCq-k 



hM 



■E 



UJ - UJk - uJ„-k + ie 



BkBq-k + 2{BkAq-k) + ^CkBq-k + 2CfcC, 



q—k 



UJ + ujk+ ujq^k - ie 



(137) 



(138) 



In the above equations we introduce the parameter A only to the calculations, 
use it as a perturbation parameter and set to one at the end of 
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If second order terms are included in the theory, they 
change the quasiparticle spectrum not only by shifting the 
quasiparticle energies but also by making them complex. The 
imaginary part that the quasiparticle energies acquire comes 
from the poles of the second order terms and it is associ- 
ated with a damping rate. The physical meaning is that when 
the energy denominator in the second order terms vanishes a 
process where a quasiparticle decays into two of lower en- 
ergy is energetically allowed. This kind of damping mech- 
anism is known as Beliaev damping and was calculated by 
Beliaev in the case of a uniform Bose superfluid ll43ll . In 
the remainder of this section we calculate the zero tempera- 
ture Beliaev damping coefficient for atoms in optical lattices 
using the tight-binding second order Beliaev approximation, 
Eqs. (fT37l-(IT38l. We follow the same ideas used by Beliaev 
to study the uniform system. 



1. Perturbative treatment 

As the starting point we assume that the net effect of sec- 
ond order plus HFB terms is to introduce small corrections to 
the Bogoliubov-de Gennes (BdG) self energies In this case 
instead of solving the equations in a self consistent way we 
can replace the BdG quasiparticle energies and amplitudes in 
the HFB and second order self-energy corrections to calculate 
the shift they introduce in the spectrum. 

The quasiparticle energies and amplitudes in the BdG ap- 
proximation are given by I1I61I : 



9.' 



(139) 



In the translationally invaiiant limit the Bogoliubov-de Gennes matrix ele- 
ments Cqq and Mq^q agree with the matrix elements calculated using the 



4(0) _ ,,(0)2 _ Eq+n'^o^U + nw'i'^ 



2 Eq + n'^o^ U — hujq^^ 



h(0) ^ ,,(0)^ ^ ^ 
'-'q ^q 



2hL0, 



(0) 



77 

C(o) - _,,(o) (0) _ "° ^ 

^1 ^ "q - „^ (0) ■ 



;(o) 



-y 



(140) 



(141) 



(142) 



(143) 



and 



g#0 



(144) 



with n the total density, n = N/M . 

As shown in the last section, the HFB approximation has 
the problem that it has a gap in the excitation spectrum and 
therefore violates Pines-Hugenholtz theorem. However, as 
shown by Beliaev It43], when second order Beliaev contribu- 
tions are included the theory becomes gapless. This can be 
seen from Eas.(fT37ll and ([138}: 

HFB-Popov \l5ttl^ 



A£q.q{0,0)-AMq-q{0,0) = -2[/m(") + — -^^ ^ 



2 /n\2\ 2 



g^7~(0) 2t/277,^ i^. ) 



(145) 



2. Beliaev damping 



If we include HFB and second order corrections, the quasi- 
particle energy shifts are given to first order in A by 
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(146) 



I 



After some algebra, Eq. \\A6\ can be written in the more con- 
venient form: 



J 



5Ea 



E 



where the matrices Bk^q-k and Bk^q-k are defined as 



,(0) 



R2 

^k,q-k 



SO) 



,(0) 



,(0) 



R2 

^k,q-k 



,(0) 



(147) 



Bk,q-k 
Bk,q-k 



'^^T^lk- 



(0) (0) 
(0) (0) 

< ^q-k - 



(0)„(0) 



(0) (0) 

< ^q-k 

(0) (0) 

")c ^q^k 



„(0)„(0) 



r 



(148) 
(149) 



If we replace Sq by hq^ /2m, the matrix elements given by 
Eqs. (I148> and ( I149> reduce to the Beliaev uniform gas matrix 
elements (see for example Ref. 1 15] and [31]). 

The damping coefficient can be obtained using the identity 
Eq. ( fT08l in Eq. ( fT47l . This yields 



7? 



Mh 



k,q-k^ 



■ OJ 



(0) 



,(0) 



(150) 

For a translational-invariant system at equilibrium, all quan- 
tities are T and R independent and depend only on the rela- 
tive coordinates r and t. Therefore, at equilibrium the scale 
separation is always valid and we can relax the condition 
Un/J < 1. In Ref. O we showed by comparison with 
solutions obtained by the exact diagonalization of the Bose- 
Hubbard Hamiltonian, that for commensurate systems in the 
parameter regime where {U/ J) < 0.5{U/ J)c, the BdG equa- 
tions give a good description of the properties of the system. 
{U / J)c ^ dn is the superfluid to Mott insulator critical ratio, 
d the dimensionality and n the density of the system. For sys- 
tems with non-commensurate fillings, where the superfluid to 
Mott insulator quantum phase transition does not take place, 
the agreement between the BdG and the exact solutions was 
shown to be significantly better for a larger parameter regime. 
Because Eq. dlSOt was found treating the second order correc- 
tions as a perturbation, its validity is restricted to the parame- 
ter regime [U / J) < 0.5{U/ J)c, where the BdG solutions are 
still a good description of the system. 

As opposed to the uniform system without the lattice, where 
for high momentum the single particle energy ( which grows 
as q^)is always dominant, in the presence of the lattice, the 
single particle excitation energies are always bounded by 4 J. 



Therefore, in the regime Un/J > 1 ^ the interaction term 
dominates for all quasimomenta and the quasiparticle ampli- 
tudes and energies can be expanded as: 



(151) 
(152) 
(153) 
where 



(0) 
Uk 


1 

2ak 




ak^ ak^ 
8 8 




1 

2ak 




ctk^ , ak^ 




~ 2n(0)C/ 


u 





a-k 



fsky/^ 



2n^°^U 



1/4 



(154) 
(155) 



In the very weakly interacting regime Urio/J < 1, the ap- 
proximations used to derive Eqs. ( I15U to jl53> are still valid 
if the quasimomentum of the excitation involved in the decay 
process is small, qai ^ y/uoU/J. 

If one substitutes Eqs. ( I151> to ( I153> for the quasiparticle 
amplitudes in Ea. (ll50> and makes use of the energy conser- 



' Notice that for large filling factors n, the parameter Un/J can be bigger 
than one but the system can be still far away from the Mott insulator critical 
point 
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vation condition, which is approximately given by 



one gets the following expression for the damping coefficient; 



7g 



97r J 



E 



eg-fc), (157) 



the one dimensional uniform Bose gas dispersion relation, it 
is also found to be zero. 

The extension of the expression for the Beliaev damping 
coefficient to higher dimensional lattice systems can be done 
straightforwardly. One just has to replace the single particle 
dispersion relation Sk in Eq |157l bv the one in the specific di- 
mension. If we assume a separable square lattice in d dimen- 
sions, with the same tunneling matrix energy J, and lattice 
constant a/, in all different directions we get: 



with Cq the dimensionless quasiparticle energies given by 
When the number of lattice sites is large. 



hull 



to a good approximation the discrete sum can be replaced by 

an integral 1/M J2k ^ «i/27r fo''^"'' dk. 

For the one dimensional system, we find that the only value 
of k at which the energy constraint is satisfied is when k = q. 
This value of k leads to a zero damping coefficient and there- 
fore in the one dimensional system the quasiparticles become 
totally stable against their decay into two of lower energy. In 
this case higher order decay processes have to be considered. 
However, the absence of Beliaev damping in one dimensional 
lattices is not a particular characteristic of the lattice disper- 
sion relation. If the damping coefficient is calculated using 



8A/'i„(o) ^ 



with the definitions £k = ^JJ2i=i^^'^^{^)' 



(0) 
k 



An analytic expression for the damping coefficient can be 
easily obtained when the excitations involved in the decay 
process have long wave number: qai ^ 1. In this parame- 
ter regime for the particular case of a three dimensional lattice 
the integral yields: 



^ J dkde sin(g) fc^ y gq£k£q-k ^ ^^^^ _ _ ai^pl _2pq COS 9^ 



327r„(0) 
9 Jaf 
32^^ 



dkk'^{q-pf 



J3 

3 



(159) 



with TO* = ti^/ {2Jaf) the effective mass. In the long wave- 
length limit, or phonon regime, the damping coefficient in the 
lattice reduces to the well known result first obtained by Beli- 
aev in the phonon regime, with the mass replaced by an effec- 
tive mass. 



Outside the phonon regime, the analytic evaluation of the 
integral is more complicated because of the energy conserva- 
tion constraint. In the uniform gas case, which has a simpler 
quasiparticle spectrum, it has been shown that there is a finite 
threshold momentum q* such that the decay of an excitation 
is impossible if q > q* |'60"]. We expect that the trigonomet- 
ric dependence on the quasimomentum of the quasiparticle 
dispersion relation in lattice-type systems makes the energy 
conservation constraint even harder to fulfill. In Ref . 1 6 1 ] the 
authors calculated the finite temperature Landau damping co- 
efficient in a one dimensional optical lattice and showed the 
disappearance of Landau damping when Uno/J > 6. 



VI. CONCLUSIONS 

In this work we continued our previous studies of the dy- 
namics of bosonic atoms confined in optical potentials. Here, 
starting from the 2PI-CTP equations of motion, derived in 
Paper 1 from the Bose-Hubbard Hamiltonian, we show how 
the complicated nonlocal, non-Markovian integro-differential 
equations can be simplified and reduced to the standard kinetic 
theory equations. Specifically, by using a two-time separa- 
tion condition, valid in dilute weakly interacting systems not 
very far away from equilibrium, we recast the full quantum 
dynamics into two coupled sets of equations: the first set of 
Boltzmann equations governing the distribution functions and 
a second set of gap equations describing the modified disper- 
sion relation. We conclude here with three remarks on some 
general features of this problem and our approach. 

First, a remark on quantum kinetic theory in discrete ver- 
sus continuous systems: Even though we work with a lattice 
gas system described by the Bose Hubbard Hamiltonian, the 
assumption that the propagators are slowly varying in the cen- 
ter of mass coordinates permits one to map the discrete tight 
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binding equations into a set of continuous differential equa- 
tions in the center of mass coordinates. For this reason the 
dynamical equations of motion we derived for discrete sys- 
tems look very similar to previous kinetic equations derived 
for continuous systems. On the other hand, to include all the 
relevant dynamical effects introduced by the lattice, we kept 
the discrete character of the tight binding Hamiltonian in the 
equations for the relative coordinates, as manifested in the gap 
equations which exhibit a dispersion relation different from 
the homogeneous Bose gas system. 

Second, the last section of this work was dedicated to a 
study of quantum equilibrium solutions. By using the quasi- 
particle approximation, we recovered from the kinetic equa- 
tions the linear HFB corrections to the self-energy plus sec- 
ond order corrections. We showed how by neglecting the 
condensate-independent second order terms in the self-energy, 
one obtains a tight-binding version of the well known Beli- 
aev equations. We used these equations to derive expressions 
for the zero temperature Beliaev damping coefficient in lat- 
tice systems in certain parameter regimes. In particular, we 
showed that for long wavelength excitations, the damping co- 
efficient in a three dimensional lattice reduces to the one cal- 
culated for a uniform Bose gas in the phonon regime, but with 
the mass replaced by the effective mass induced by the lattice. 

A final remark on the purpose of this work. It is not meant 



to be a mere academic exercise in our demonstration of how 
Boltzmann like equations are obtained from the effective ac- 
tion and equilibrium solutions can be obtained from the full 
quantal solutions. In making explicit the simplifying assump- 
tions en route starting from first principles, it allow us to re- 
alize the limitations and the applicability of a kinetic theory 
formulation for describing the quantum dynamics of many- 
body lattice systems. It serves to identify the range of validity 
and the parameter regimes where the underlying assumptions 
leading to these simplified kinetic equations can become un- 
reliable. We view this effort as having both theoretical and 
practical significance in seeking a proper description of such 
systems and better understanding of its behavior - theoretical 
in scrutinizing the practicing kinetic theories in existence, and 
practical in providing the correct paramters for comparison 
with experiments. 
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